Finite temperature lattice properties of graphene beyond the quasiharmonic 

approximation 
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The thermal and mechanical stability of graphene is important for many potential applications in 
nanotechnology. We calculate the temperature dependence of lattice parameter, elastic properties 
and heat capacity by means of atomistic Monte Carlo simulations that allow to go beyond the quasi- 
harmonic approximation. We predict an unusual, non-monotonic, behavior of the lattice parameter 
with minimum at T « 900 K and of the shear modulus with maximum at the same temperature. 
The Poisson ratio in graphene is found to be small m 0.1 in a broad temperature interval. 



PACS numbers: 81.05.Uw, 62.20. 
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Understanding the structural and thermal properties 
of two dimensional (2D) systems is one of the challenging 
problems in modern statistical physics [l| . Traditionally, 
it was discussed mainly in the context of biological mem- 
branes and soft condensed matter. The complexity of 
these systems hindered any truly microscopic approach 
based on a realistic description of interatomic interac- 
tions. The discovery of graphene [2|, the first truly 2D 
crystal made of just one layer of carbon atoms, provides 
a model system for which an atomistic description be- 
comes possible. The interest for graphene has been trig- 
gered by its exceptional electronic properties (for review 
see [1, 0, Hi) but the experimental observation of ripples 
in freely hanged graphene Q has initiated a theoretical 
interest also in the structural properties of this mate- 
rial 0,1- Ripples or bending fluctuations have been pro- 
posed as one of the dominant scattering mechanisms that 
determine the electron mobility in graphene 0. More- 
over, the structural state influences the mechanical prop- 
erties that are important in themselves for numerous po- 
tential applications of graphene [TH EH E3] • 

Two dimensional crystals are expected to be strongly 
anharmonic due to an intrinsic bending instability cou- 
pled to in-plane stretching modes. This coupling is cru- 
cial to prevent crumpling of the crystal and stabilize 
the flat phase These expectations have been con- 
firmed by atomistic simulations for graphene showing 
very strong bond length fluctuations already at room 
temperature H. Beside the relevance for 2D systems, 
anharmonicity [131 ] is of general importance in condensed 
matter in relation to structural phase transitions (T3.[l5|. 
soft modes in ferroelectrics [l6| , melting [l7| and related 
phenomena. Usually anharmonicity in crystals is weak 
enough and thus can be well described in the frame- 
work of perturbation theory [ID, [ll| 0, [2(|. However, 
this might be not the case for strongly anharmonic sys- 
tems, like graphene. Atomistic simulations offer the pos- 
sibility to study anharmonic effects for a specific ma- 
terial without need of perturbative schemes. For car- 
bon a very accurate description of energetic and ther- 
modynamic properties of different allotropes including 
graphene 0, [2l[ is provided by the empirical bond order 
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FIG. 1: (color online) Temperature dependence of the lattice 
parameter a (solid blue line) and nearest neighbor distance 
Run (dashed red line). The scales of left (a) and right (Rnn) 
y-axes are related to each other by v3- At T = 0, a — 
2.4595 • 10~ 10 m. 



potential LCBOPII 22]. Here we present the tempera- 
ture dependence of thermodynamical and elastic proper- 
ties of graphene, calculated by means of atomistic Monte 
Carlo (MC) simulations based on LCBOPII. 

We perform MC simulations at finite temperature T 
with periodic boundary conditions for a sample of N — 
8640 atoms with equilibrium size at zero temperature 
of 147.57 A in the x direction and 153.36 A in the y 
direction. We equilibrate the sample in the NPT en- 
semble at pressure P = for at least 2 • 10 5 MC steps 
(1 MC step corresponds to N attempts to a coordinate 
change) which we found to be enough for convergence 
of total energy and sample size. Further 10 5 MC steps 
are used to evaluate the average lattice parameter a and 
average nearest neighbor distance R nn and radial distri- 
bution function g(R). 

Figure [1] shows that a and R nn decrease with increas- 
ing temperature up to about 900 K, yielding a negative 
thermal expansion coefficient a = (— 4.8± 1.0) • 10~ 6 K _1 
in the range 0-300 K. As noted in Ref. [23[ this anomaly 
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FIG. 2: a) (color online) Nearest neighbor radial distribution 
function g(R nn ) for the N = 8640 sample at 300 K, 900 K 
and 2300 K. The vertical lines indicate the length of double 
(1.34- 10" 10 m), conjugated (1.42- 10" 10 m) and single (1.54- 
10 -10 m) bonds, b) Standard deviation a(R„ n ) (red circles) 
and the best fit to yT in the temperature range up to 500 K 
(solid blue line). 



FIG. 3: a) (color online) 2D elastic moduli of graphene as 
of function of temperature: adiabatic bulk modulus i>A (solid 
blue line with circles), isothermal bulk modulus 6t (dashed 
blue line with squares) and shear modulus fj, (solid red line 
with diamonds), b) Adiabatic Young's modulus Ya (solid 
blue line with circles) and isothermal Yt (dashed red line 
with squares). 



is due to a low- lying bending phonon branch [24| ■ Our re- 
sults are in a gree ment up to 500 K with those of Mounet 
and Marzari [23| who used the quasi-harmonic approxi- 
mation with phonon frequencies and Griineisen parame- 
ters calculated by the density functional approach. How- 
ever, at higher temperatures our results are qualita- 
tively different, since in Ref. (23j a remains negative in 
the whole studied temperature interval up to 2200 K, 
whereas we find that it changes sign and becomes pos- 
itive at T w 900 K. This discrepancy with the quasi- 
harmonic theory, which in general works reasonably well 
for three-dimensional crystals, is one of the evidences of 
strong anharmonicity in graphene. 

The deviations from harmonic behavior can be char- 
acterized by examining the radial distribution function 
g(R) around the first neighbor distance R nn — 1.42 • 
10 -10 m. In Fig. [2Ja) we present g(R) and the related 
standard deviation a{R nn ) shown in Fig. HJb). In the 
harmonic approximation R nn would have a Gaussian dis- 
tribution yilding a(R nn ) oc y/T. Deviations from square 
root behavior can be observed above 900 K, achieving 
10 % at 2000 K. 

The Lindemann criterion has been shown to apply also 
in 2D, giving a(R nn ) ~ 0.23i? rm at melting 25]. We 
found a(R nn )/R nn = 0.056 at T = 2300 K, indicating 
that we are significantly below melting point. Moreover, 
conventional theory of two-dimensional melting relates it 
to the formation of topological defects [26| . In our simu- 
lations we have not seen any sign of premelting anomalies 
(formation of vacancies, topological defects etc.) up to 
3500 K 0. 

The strong anharmonic behavior of graphene leads also 
to unusual temperature dependence of the elastic moduli. 
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FIG. 4: Adiabatic va and isothermal vt Poisson ratio of 
graphene as function of temperature. 

The 2D bulk modulus b is defined by 

E is - 2bui, (I) 

where Ei S is the elastic energy per unit area under an 
isotropic deformation u yy = u xx = u, s , u xy = 0. 

For uniaxial deformations u xx (u yy = u xy = 0) the 
elastic energy is 

Eum = \(b + v)u 2 xxl (2) 

where /i is the 2D shear modulus. 

Isothermal moduli are also expressed as in Eq. (fTJ) and 
Eq. ©, with replacement of the energy E by the free 
energy F = —TlnZ where Z is the partition function. 
Although it is impossible in MC to calculate F directly, 
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we will use the fact that adiabatic and isothermal shear 
moduli /i coincide [27| and that the Poisson ratio defined 
below can be calculated directly to derive the isothermal 
bulk modulus bx- The Young modulus Y and Poisson 
ratio v are defined in terms of b and /j, as [28| 



Y 



Abfi 
b + fi' 

b — fi 
6 + At" 



(3) 
(4) 



The Poisson ratio can also be defined as the ratio be- 
tween the axial e ax i a i and transverse e tra ns strain as 



(5) 



TABLE I: Adiabatic bulk (6a), shear (/j.) and isothermal bulk 
(for) moduli, and isothermal Poisson ratio (vt)- 



T, K b A (eV • A" 2 ) fj, (eV • A~ 2 ) 
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0.16i0.03 
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9.80i0.15 
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300 


12.36±0.04 


9.95i0.17 


12.52il.41 


0.12i0.05 


500 


12.22±0.05 


10.16i0.20 


12.24il.66 


0.09i0.06 


700 


12.09±0.05 


10.27i0.17 


12.93i2.13 


0.12i0.08 


900 


11.94T0.04 


10.25i0.18 


11.29i2.20 


0.09i0.09 


1100 


11.85±0.06 


10.21i0.22 


11.31i2.57 


0.05i0.11 


1300 


11.70±0.04 


10.07i0.21 


12.05i3.00 


0.09i0.12 


1500 


11.57T0.04 


9.94i0.18 


11.63i3.10 


0.08i0.13 


1700 


11.44T0.04 


9.75i0.24 


8.44i3.20 


-0.07i0.18 


2000 


11.31i0.06 


9.52i0.22 






2100 


11.23i0.05 


9.46i0.26 


8.26i3.58 


-0.07i0.21 



The latter definition provides a way to calculate the 
isothermal vt so that Eq. ([¥]) with \ia = (J-t = fJ- 
yields br- 

Adiabatic bulk and shear moduli 6 a and /j, have been 
calculated using the following procedure. We equilibrate 
the sample as described before. Afterwords, 20 config- 
urations separated by 5000 MC steps were stored and 
subjected to either isotropic or uniaxial deformation in 
steps of 0.01 % without letting the sample relax. For each 
sample, the variation of the elastic energy with deforma- 
tion was then fitted to Eq. (fT]) and Eq. ^ over 21 points 
around the undistorted configuration. The averages of 
the calculated bA and [i for the 20 samples are given in 
Table I and shown in Fig. [^together with the derived Ya- 
We find that the temperature dependence of [i is anoma- 
lous. While in general all elastic moduli decrease as a 
function of temperature due to weakening of interatomic 
interactions with temperature, in graphene /i grows with 
increasing temperature up to T ~ 700 — 900 K which is 
the same temperature where the thermal expansion be- 
havior (Fig. [1]) becomes normal. The Young modulus 
Y follows the same anomalous temperature dependence 
as /i. 

We find that the behavior of the elastic energy as a 
function of deformation u is parabolic in a wider range 
of deformations, up to about 0.2 %. For larger deforma- 
tions, the elastic energies follow a cubic dependence on 
the deformation at least up to u = 3 %. At this value 
the ratio of the cubic term to the quadratic one in the 
elastic energy is about 0.12. Up to 10 % deformation and 
up to 2200 K, deformations are reversible, and no defect 
(vacancy and Stone- Wales [2!| or dislocations [3(|) are 
found. This is not surprising in view of the very high 
cohesive energy (7.6 eV/atom in graphite [22] ) of carbon 
and defect formation energy in graphene [29j |. To the 
best of our knowledge, there are no experimental data 
on defect formation under strain in this range of temper- 
atures. 

Next, the isothermal Poisson ratio vt has been cal- 
culated using the following procedure. We take the 



graphene sample equilibrated as described before at a 
given temperature. The sample is then stretched of 1 % 
in the x and y directions separately and re-equilibrated 
again for at least 5- 10 4 MC steps. After re-equilibration, 
the sample size in the x and y directions have been aver- 
aged for at least 5 • 10 4 MC steps and the corresponding 
strain e x and e y have been calculated yielding the Pois- 
son ratio in each direction through Eq. ([5]) . The Poisson 
ratios in the x and y directions are very close and we 
take their average as vt- The calculated adiabatic and 
isothermal Poisson ratios va and vt, shown in Fig. 0] 
and Table I, are very small and coincide within the er- 
ror in the whole studied temperature range. However 
at high temperature, we find that vt can become neg- 
ative. Materials with negative Poisson ratio are called 
auxetic and, in general, this property is related to very 
unusual crystalline structures. Membranes, on the other 
hand, may display this behavior due to entropy. In fact, 
an expansion in the unstretched direction contrasts the 
reduction of phase space due to the decrease of height 
fluctuations due to stretching. Furthermore, the small- 
ness of v implies that the Lame' constant A = b — fi 
is small in comparison with //. Therefore for a generic 
deformation described by a tensor u, the elastic energy 
E e i = Ilu\a + (l/2)A(Tru 2 ) 1] for graphene can be ap- 
proximated as E e i w 1-iiUij. 

Once vt is known we can calculate &t from Eq. (|U) 
and Yt from Eq. ([3]). The calculated &t and Yr are 
presented in Table I and compared to the adiabatic values 
in Fig. El At T = 300 K, we find Ya — 353 i 4 N • ur 1 
and Yt — 355 ± 21 N • m _1 in good agreement with the 
experimental value 340 i 50 N • m -1 [ill ]. 

Another important anharmonic effect is the tempera- 
ture dependence of the molar heat capacity at constant 
volume C V (T) = 3R(1 + T/E ) (see Fig. EJ), where R is 
the gas constant and E is a typical energy of interatomic 
interactions (l3l . [20| . The low temperature behavior was 
calculated in the harmonic approximation in [23]. Our 
approach is classical and therefore can be used to calcu- 
late Cy only at high temperatures. On the other hand, 
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FIG. 5: Temperature dependence of the molar heat capacity 
at constant volume Cv (solid line) and fit to Cv(T) — 3i?(l + 
T/E ) (dashed line). 



our approach does not use the harmonic approximation, 
yielding information about phonon-phonon interaction 
effects. Our calculations show that the linear temper- 
ature dependence of Cy becomes noticeable for T > 800 
K with Eg = 1.3 eV. Contrary to alkali metals where Eq 



is of the order of the vacancy formation energy [l8[, for 
graphene, due to anharmonicity, Eq is about 1/5 of the 
defect formation energy. 

In summary we have presented the temperature depen- 
dence of lattice parameter, elastic moduli and high tem- 
perature heat capacity of graphene calculated by Monte 
Carlo simulations based on the LCBOPII empirical po- 
tential 22] for a crystallite of about 15 x 15 nm 2 . In the 
studied range of temperatures, up to 2200 K, and for de- 
formation as large as 10% we have not seen any sign of 
defect formation. Indeed the very high energy for defect 
formation in graphene makes this material exceptionally 
strong, as also found experimentally [HI, [HJ. We find 
that graphene is strongly anharmonic due to soft bend- 
ing modes yielding strong out of plane fluctuations. We 
find that, up to 900 K, graphene is anomalous since its 
lattice parameter decreases and shear modulus increases 
with increasing temperature going over to normal be- 
havior at higher temperatures. It would be interesting to 
check these predictions experimentally. 
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Wetenschappelijk Onderzoek (NWO)'. 
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